*** Code for Remotely Incorrect? *****
******************************
***Misclassification Models***
***********Ejido FEs**********
******************************

clear all
set more off
cap log cl
cap prog drop mclogit
cap prog drop mclogit_lf
cap prog drop mclogit_p
cap prog drop mcre
cap prog drop mcre_p
set seed 2356935

local 	sel_user		2		// 1 - Dan, 2 - Jen
	
	if 		`sel_user' 		== 1 { 
		global 	user		"C:\Users\03638881\Dropbox\Satellite data measurement\"
	}	
	
	if 		`sel_user' 		== 2 {
		global user 		"C:\Users\alixgarj\Dropbox\Satellite data measurement\"
	}
	

	
	global	data			"${user}Data\"
	global 	outputs			"${user}Latex\tables\"
	global 	graphs	 		"${user}Latex\graphs"
	global 	dofiles			"${user}Do files\"

adopath + "${user}Do Files\"


log using "${outputs}MeasurementErrori2iDataJAG2", replace
estimates clear

**Spec 1: Full Sample, Ejido FEs**
use "${data}PESi2iPanelDataShortAugust2021.dta", clear
qui {
	loc spec = 11
	keep if inejido==1
	loc fe "EjidoCode"
	loc cluster "EjidoCode"
	sort APEJNCNID EjidoCode yr
	xtset `fe'

	*define vars to be used
	loc y "binary_defor"

	xi i.yr
	loc x0 "_Iyr_2002 _Iyr_2003 _Iyr_2004 _Iyr_2005 _Iyr_2006 _Iyr_2007 _Iyr_2008 _Iyr_2009 _Iyr_2010 _Iyr_2011 _Iyr_2012 _Iyr_2013 _Iyr_2014"
	loc x1 "beneL1 everbene elev_met_mean slope_deg_mean dist_any_road_mt dist_5000_city_km TPApAreaHa pctTC2000 per_maj_ind"

	*rescale vars*
	replace TPApAreaHa=TPApAreaHa/100
	replace dist_any_road_mt=dist_any_road_mt/1000
	replace dist_5000_city_km=dist_5000_city_km/100
	replace slope_deg_mean=slope_deg_mean/100
	replace elev_met_mean=elev_met_mean/100

	g z1=path_row_counts_2000_2015_L7
	g z2=path_row_counts_2000_2015_L7*slope_deg_mean 
	g z3=path_row_counts_2000_2015_L7*TPApAreaHa
	lab var path_row_counts_2000_2015_L7 "Cloud-free scenes, L7"
	lab var z2 "Cloud-free scenes, L7 x Avg Slope"
	lab var z3 "Cloud-free scenes, L7 x Area"
	lab var dist_5000_city_km "Distance to city with $>$ 5,000 people"

	loc z1 "z1"
	loc z2 "z2 z3"

	*create means for CRE models*
	loc mx1 " "
	loc mz1 " "
	loc mz2 " "
	foreach x in `x1' {
		bys `fe': egen m`x'=mean(`x')
		loc mx1 = "`mx1' m`x'"
	}
	foreach x in `z1' {
		bys `fe': egen m`x'=mean(`x')
		loc mz1 = "`mz1' m`x'"
	}
	foreach x in `z2' {
		bys `fe': egen m`x'=mean(`x')
		loc mz2 = "`mz2' m`x'"
	}
	
	loc x1 "elev_met_mean slope_deg_mean dist_any_road_mt dist_5000_city_km TPApAreaHa pctTC2000 per_maj_ind"

}

*estimators: misclassification*
*mc logit**
loc i=1

logit `y' `x0' i.beneL1 i.everbene `x1' `mx1', vce(cl `cluster') nolog
loc b=_b[_cons]
loc fromp "`y':_cons=`b'"
foreach x in 1.beneL1 1.everbene `x1' {
	loc b=_b[`x'] 
	loc fromp "`fromp' `y':`x'=`b'"
}
loc fromp "`fromp' a0:_cons=-2.5 a0:z1=-2.5 a1:_cons=-.4 a1:z1=.05"
mclogit `y' i.beneL1 i.everbene `x1', vce(cl `cluster') mc(`z1') diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in 1.beneL1 1.everbene `x1' {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}
foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1
	loc fromp "`fromp' a`q':_cons=`b1' a`q':z1=`b2'"		
}
mclogit `y' `x0' i.beneL1 i.everbene `x1', vce(cl `cluster') mc(`z1') diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in `x0' 1.beneL1 1.everbene `x1' {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}
foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1
	loc fromp "`fromp' a`q':_cons=`b1' a`q':z1=`b2'"		
}

mclogit `y' `x0' i.beneL1 i.everbene `x1' `mx1', vce(cl `cluster') mc(`z1') diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
mat b1=e(b),0
loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in `x0' 1.beneL1 `x1' `mx1' {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}
foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1		
}
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
foreach x in `z1' {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
eststo col`i': margins, dydx(beneL1 everbene `x1') predict(pr) post
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1

loc fromp "`fromp' a0:_cons=-1.5 a0:z1=-.2 a0:z2=.15 a0:z3=.002 a1:_cons=-.5 a1:z1=-.3 a1:z2=.3 a1:z3=.005"		
mclogit `y' `x0' i.beneL1 i.everbene `x1' `mx1', vce(cl `cluster') mc(`z1' `z2') diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
mat b2=e(b),0
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
g G0l_`spec'=[a0]_cons
g G1l_`spec'=[a1]_cons
foreach x in `z1' `z2' {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
	replace G0l_`spec'=G0l_`spec'+`x'*[a0]`x'
	replace G1l_`spec'=G1l_`spec'+`x'*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
replace G0l_`spec'=normal(G0l_`spec')
replace G1l_`spec'=normal(G1l_`spec')

cap drop xb
predict xb, xb
g mcl_`spec'_beneL1 = logistic(xb+_b[1.beneL1]*(1-beneL1)) - logistic(xb-_b[1.beneL1]*beneL1) if e(sample)
g mcl_`spec'_everbene = logistic(xb+_b[1.everbene]*(1-everbene)) - logistic(xb-_b[1.everbene]*everbene) if e(sample)				
foreach x in `x1' {
	g mcl_`spec'_`x' = logisticden(xb)*[`y']`x' if e(sample)
}

eststo col`i': margins, dydx(beneL1 everbene `x1') predict(pr) post 
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1	

*estimators: misclassification*
*mc scobit*


loc logL = -999999999

forval a=  0.2(0.05)0.95 {

	loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	constraint 2 [ln_a2]_cons = `aa'
	constraint 3 [a0]_cons = -10
	constraint 4 [a0]z1 = 0
	constraint 5 [a0]z2 = 0
	constraint 6 [a0]z3 = 0
	constraint 7 [a0]z4 = 0
	
cap	mcre `y' `x0' i.beneL1 i.everbene `x1' `mx1', vce(cl `cluster') mc(`z1') diff from(b1, copy) constraint(2 3 4) tech(nr 5 dfp 30 bhhh 10) nolog
	if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				loc ll = e(ll)
				loc a0=[a0]_cons
				loc a1=[a1]_cons
				foreach x in `z1' {
					su `x'
					loc a0=`a0'+r(mean)*[a0]`x'
					loc a1=`a1'+r(mean)*[a1]`x'
					
					}
					loc a0=normal(`a0')
					loc a1=normal(`a1')
					eststo col`i': margins, dydx(beneL1 everbene `x1') predict(pr) post
					eststo col`i', add(G0 `a0')
					eststo col`i', add(G1 `a1')
					eststo col`i', add(alpha `a')
					eststo col`i', add(logL `ll')
			}
	}
	
}
	loc i=`i'+1


loc logL = -999999999	
forval a=  0.2(0.05)0.95 {
loc h=round(`a'*100,5)
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	constraint 2 [ln_a2]_cons = `aa'
	constraint 3 [a0]_cons = -10
	constraint 4 [a0]z1 = 0
	constraint 5 [a0]z2 = 0
	constraint 6 [a0]z3 = 0
	constraint 7 [a0]z4 = 0
	cap mcre `y' `x0' i.beneL1 i.everbene `x1' `mx1', vce(cl `cluster') mc(`z1' `z2') diff from(b2, copy) constraint(2 3 4 5 6) tech(nr 5 dfp 30 bhhh 10) nolog
	if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				loc ll = e(ll)
				loc a0=[a0]_cons
				loc a1=[a1]_cons
				g G0s`h'_`spec'=[a0]_cons
				g G1s`h'_`spec'=[a1]_cons
				foreach x in `z1' `z2' {
					su `x'
					loc a0=`a0'+r(mean)*[a0]`x'
					loc a1=`a1'+r(mean)*[a1]`x'
					replace G0s`h'_`spec'=G0s`h'_`spec'+`x'*[a0]`x'
					replace G1s`h'_`spec'=G1s`h'_`spec'+`x'*[a1]`x'
					}
				loc a0=normal(`a0')
				loc a1=normal(`a1')
				replace G0s`h'_`spec'=normal(G0s`h'_`spec')
				replace G1s`h'_`spec'=normal(G1s`h'_`spec')

				*cap drop xb
				*predict xb, xb
				loc alp = exp([ln_a2]_cons)
				g mcs_`h'_`spec'_beneL1 = 1 - (1/(1+exp(xb+_b[1.beneL1]*(1-beneL1)))^`alp') - ( 1 - (1/(1+exp(xb-_b[1.beneL1]*beneL1))^`alp') ) if e(sample)
				g mcs_`h'_`spec'_everbene = 1 - (1/(1+exp(xb+_b[1.everbene]*(1-everbene)))^`alp') - ( 1 - (1/(1+exp(xb-_b[1.everbene]*everbene))^`alp') ) if e(sample)				
				foreach x in `x1' {
					g mcs_`h'_`spec'_`x' = [`y']`x'*exp(xb)*`alp'*((1+exp(xb))^((-1)*`alp' - 1)) if e(sample)
					}

	eststo col`i': margins, dydx(beneL1 everbene `x1') predict(pr) post
	eststo col`i', add(G0 `a0')
	eststo col`i', add(G1 `a1')
	eststo col`i', add(alpha `a')
	eststo col`i', add(logL `ll')
	
	
		}
	}
}

loc i=`i'+1	



	*Partial Observability Probit
g binary_defor1=binary_defor
probit binary_defor `x0' i.beneL1 i.everbene `x1' `mx1' if `z1'!=.
loc b=_b[_cons]
loc fromp "binary_defor:_cons=`b'"
foreach x in `x0' 1.beneL1 1.everbene `x1' `mx1' {
	loc b=_b[`x']
	loc fromp "`fromp' binary_defor:`x'=`b'"
}
probit binary_defor1 `z1' `z2'
loc b=_b[_cons]
loc fromp "`fromp' binary_defor1:_cons=`b'"
foreach x in `z1' `z2' {
	loc b=_b[`x']
	loc fromp "`fromp' binary_defor1:`x'=`b'"
}
	
biprobit (binary_defor `x0' i.beneL1 i.everbene `x1' `mx1') (binary_defor1 `z1' `z2'), partial from(`fromp') diff tech(nr 25 dfp 25) nolog 
loc ll=e(ll)
loc a1=[binary_defor1]_cons
g G1po_`spec'=[binary_defor1]_cons
foreach x in `z1' `z2' {
	su `x'
	loc a1=`a1'+r(mean)*[binary_defor1]`x'
	replace G1po_`spec'=G1po_`spec'+`x'*[binary_defor1]`x'
}
loc a1=1-normal(`a1')
replace G1po_`spec'=1-normal(G1po_`spec')

cap drop xb
predict xb, xb1
g mcpo_`spec'_beneL1 = normal(xb+_b[1.beneL1]*(1-beneL1)) - normal(xb-_b[1.beneL1]*beneL1) if e(sample)
g mcpo_`spec'_everbene = normal(xb+_b[1.everbene]*(1-everbene)) - normal(xb-_b[1.everbene]*everbene) if e(sample)
foreach x in `x1' {
	g mcpo_`spec'_`x' = [binary_defor]`x'*normalden(xb) if e(sample)
}

eststo col`i': margins, dydx(beneL1 everbene `x1') pr(pmarg1) post	
eststo col`i', add(G0 0)
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1	

esttab col*  ///
	using "${outputs}i2iresults11po.tex", ///
	se nonotes  style(tex)  b(%12.3f) se(%12.3f)  noobs ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" "" "" "") ///
	 keep(1.beneL1 1.everbene `x1') nonumbers  replace  fragment  ///
	 scalar(G0 G1 alpha logL) prehead({ \begin{tabular}{lccccc}    ///
	 \hline \hline ///
	 & \multicolumn{2}{c}{MC Logit} & \multicolumn{2}{c}{MC Scobit} & \multicolumn{1}{c}{Partial Observability} \\  ///
	 & (1) & (2) & (3) & (4) & (5) \\ ) ///
	 prefoot( \\ ) ///
		postfoot(Full MC covars & no & yes & no & yes & yes \\ \hline \hline \end{tabular} }  ///
		\begin{tablenotes}[para,flushleft] \footnotesize{Column headers indicate estimator. ///
		Standard errors are clustered at the municipality level. ///
		Marginal effects evaluated at sample means are displayed. $G_{0}$ and $G_{1}$ are the probability of a false positive and negative, respectively, ///
		evaluated at sample means. Column 1 allows the misclassification rates to depend ///
		on the number of L7 cloud-free scenes. Column 2 allows the misclassification rates to depend on the number of L7 cloud-free scenes and its ///
		interaction with average slope and polygon area. Time fixed effects included in all models. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )

estimates clear

keep `y' beneL1 everbene `x1' `z1' `z2' yr muni_cve i2itpapid EjidoCode APEJNCNID G0l* G1l* mcl* G0s* G1s* mcs* G1po* mcpo*
save "${data}PESi2iPanelData_results_mc_po_ej.dta", replace



	estimates clear



log cl